
#This project is statistical analysis for a 3X1 factorial experiment-1 and 3x2 factorial experiment-2. For experiment-1 contain four scenarios for each of three conditions.
# The conditions involve AI, Human and Human aided by AI. each of these are dummy coded # to run the regression model.
# TRUST, CREDEBILITY, AGREEABLENESS, AND FAIRNESS are the dependent variables for experiment-1. 
# TRUST, CREDEBILITY, QUALITY, AND BIAS are the dependent variables for experiment-2 as perceived by #participants. conditions is independent variable.

# Step-0: pre-process data to remove participants who does not qualify
# Step-1: regression analysis for experiment-1 and 2
# Step-2: MLM models experiment-1.
# Step-2: MLM models experiment-2.

########################################################################
#                                                                      #
#                     0 Data cleaning process                          #
#                                                                      #
########################################################################

# Loading relevant libararies
library(MASS)
library(dplyr)
library(scales)
library(performance)
library(reshape2)
library(nlme)
library(car)
require(ggplot2)
require(ggpubr)

# load the data (data/directory name go below.)
ex.data=read.csv("xxxxxxxxxxxxxxxxxxxxxx")

# filter participants who failed attention check and not active online
dat_com = ex.data %>% filter((SOCIAL_1!=1 | SOCIAL_2!=1 | SOCIAL_3 !=1 | SOCIAL_4!=1 ) & CHECK==2)

# Find 48% of the median time participants took to complete the experiment
dat_com$Duration..in.seconds. <- as.numeric(dat_com$Duration..in.seconds.)
class(dat_com$Duration..in.seconds.)
dat_com_str<- dat_com
median(as.numeric(ex.data$Duration..in.seconds.))
quantile(ex.data$Duration..in.seconds., na.rm = T)

#excluding-speeders i.e. participants who took less than 48% of the median time
dat_com_str_sp<-dat_com_str %>%
  filter(Duration..in.seconds. > 295.44)
nrow(dat_com_str_sp)

# excluding duplicates
final_data_Spain <- dat_com_str_sp %>% distinct(pid, .keep_all = T)
nrow(final_data_Spain)

# excluding participants who did not complete the experiment
final_data_Spain <- final_data_Spain %>% filter(!is.na(final_data_Spain$BEN_CRED))
nrow(final_data_Spain)


##################################################################################################
#                                                                                                #
#                   1 Regression analysis for experiment-1                                       #
#                                                                                                #
##################################################################################################
########################################################################
#                                                                      #
#               1.1 Creating variables of interest                     #
#                                                                      #
########################################################################

#### creating dummy variables (IVs) for AI, Humnan, and Mix conditions
final_data_Spain$AI.dummy <- ifelse(final_data_Spain$condition %in% c(2,3),0,1)
final_data_Spain$Human.dummy <- ifelse(final_data_Spain$condition %in% c(1,3),0,1)
final_data_Spain$AI_Human.dummy <- ifelse(final_data_Spain$condition %in% c(1,2),0,1)

#### Creating aggregated dependent variables (trust, fair, bias, and agree) across scenarios
final_data_Spain$sp1mean.trust <- (final_data_Spain$TRUST+final_data_Spain$CIVIL_TRUST+ final_data_Spain$FACTS_TRUST+ final_data_Spain$DEC_TRUST)/4
final_data_Spain$sp1mean.fair <- (final_data_Spain$FAIR+final_data_Spain$CIVIL_FAIR+ final_data_Spain$FACTS_FAIR+ final_data_Spain$DEC_FAIR)/4
final_data_Spain$sp1mean.biased <- (final_data_Spain$BIASED+final_data_Spain$CIVIL_BIASED+ final_data_Spain$FACTS_BIASED+ final_data_Spain$DEC_BIASED)/4
final_data_Spain$sp1mean.agree <- (final_data_Spain$AGREE+final_data_Spain$CIVIL_AGREE+ final_data_Spain$FACTS_AGREE+ final_data_Spain$DEC_AGREE)/4

### Creating aggregated dependent variables for each scenarios

#scenario1
final_data_Spain$agg_s1 <- (final_data_Spain$TRUST + final_data_Spain$FAIR + final_data_Spain$AGREE + final_data_Spain$BIASED)/4 

#scenario2
final_data_Spain$agg_s2 <- (final_data_Spain$CIVIL_TRUST + final_data_Spain$CIVIL_FAIR + final_data_Spain$CIVIL_AGREE + final_data_Spain$CIVIL_BIASED)/4 

#scenario3
final_data_Spain$agg_s3 <- (final_data_Spain$FACTS_TRUST + final_data_Spain$FACTS_FAIR + final_data_Spain$FACTS_AGREE + final_data_Spain$FACTS_BIASED)/4  

#scenario4
final_data_Spain$agg_s4 <- (final_data_Spain$DEC_TRUST + final_data_Spain$DEC_FAIR + final_data_Spain$DEC_AGREE + final_data_Spain$DEC_BIASED)/4


########################################################################
#                                                                      #
#    1.2 Randomization check across different control variables        #
#                                                                      #
########################################################################

#### randomization check is performed to determine that participants were randomly distributed across the conditions
#summary(aov(scale(condition)~scale(RACE), data = subset(final_data, !is.na(RACE)) ), na.rm=T)
summary(aov(scale(condition)~scale(INCOME) , data = subset(final_data_Spain,!is.na(INCOME))), na.rm=T)
summary(aov(scale(condition)~scale(GENDER) , data = subset(final_data_Spain,!is.na(GENDER))), na.rm=T)
summary(aov(scale(condition)~scale(EDU), data = subset(final_data_Spain, !is.na(EDU))), na.rm=T)
summary(aov(scale(condition)~scale(AGE_CAT), data = subset(final_data_Spain, !is.na(AGE_CAT))), na.rm=T)
#summary(aov(scale(condition)~scale(PARTY7), data = subset(final_data, !is.na(PARTY7))), na.rm=T)
summary(aov(scale(condition)~scale(IDEO_1), data = subset(final_data_Spain, !is.na(IDEO_1))), na.rm=T)
summary(aov(scale(condition)~scale(FOLLOW), data = subset(final_data_Spain, !is.na(FOLLOW))), na.rm=T)
summary(aov(scale(condition)~scale(INCOME), data = subset(final_data_Spain, !is.na(FOLLOW))), na.rm=T)


########################################################################
#                                                                      #
# 1.3 Regression analysis for study1 (aggregated dependent variables)  #
#                                                                      #
########################################################################

summary(aov((sp1mean.trust) ~ (condition), data = final_data_Spain), na.rm=T)
summary(aov((sp1mean.fair) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((sp1mean.biased) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((sp1mean.agree) ~ (condition), data = final_data_Spain), na.rm=T)

summary(lm((sp1mean.trust) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((sp1mean.fair) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((sp1mean.biased) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((sp1mean.agree) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

summary(lm((sp1mean.trust) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((sp1mean.fair) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((sp1mean.biased) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((sp1mean.agree) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

########################################################################
#                                                                      #
#   1.4 Regression analysis for study1 (Scenariowise for each dv)      #
#                                                                      #
########################################################################

################################################### Scenario-1 ######################################################
summary(aov((TRUST) ~ (condition), data = final_data_Spain), na.rm=T)
summary(aov((FAIR) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((BIASED) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((AGREE) ~ (condition), data = final_data_Spain), na.rm=T)

summary(lm((TRUST) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FAIR) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((BIASED) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((AGREE) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

summary(lm((TRUST) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FAIR) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((BIASED) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((AGREE) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

################################################### Scenario-2 ######################################################
summary(aov((CIVIL_TRUST) ~ (condition), data = final_data_Spain), na.rm=T)
summary(aov((CIVIL_FAIR) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((CIVIL_BIASED) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((CIVIL_AGREE) ~ (condition), data = final_data_Spain), na.rm=T)

summary(lm((CIVIL_TRUST) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)
summary(lm((CIVIL_FAIR) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)
summary(lm((CIVIL_BIASED) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((CIVIL_AGREE) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

summary(lm((CIVIL_TRUST) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)
summary(lm((CIVIL_FAIR) ~ (Human.dummy) + (AI_Human.dummy) , data =final_data_Spain), na.rm=T)
summary(lm((CIVIL_BIASED) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)
summary(lm((CIVIL_AGREE) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

################################################### Scenario-3 ######################################################
summary(aov((FACTS_TRUST) ~ (condition), data = final_data_Spain), na.rm=T)
summary(aov((FACTS_FAIR) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((FACTS_BIASED) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((FACTS_AGREE) ~ (condition), data = final_data_Spain), na.rm=T)

summary(lm((FACTS_TRUST) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FACTS_FAIR) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FACTS_BIASED) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FACTS_AGREE) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

summary(lm((FACTS_TRUST) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FACTS_FAIR) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FACTS_BIASED) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((FACTS_AGREE) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

################################################### Scenari-4 ######################################################
summary(aov((DEC_TRUST) ~ (condition), data = final_data_Spain), na.rm=T)
summary(aov((DEC_FAIR) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((DEC_BIASED) ~ (condition) , data = final_data_Spain), na.rm=T)
summary(aov((DEC_AGREE) ~ (condition), data = final_data_Spain), na.rm=T)

summary(lm((DEC_TRUST) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((DEC_FAIR) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((DEC_BIASED) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((DEC_AGREE) ~ (AI.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

summary(lm((DEC_TRUST) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((DEC_FAIR) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((DEC_BIASED) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((DEC_AGREE) ~ (Human.dummy) + (AI_Human.dummy), data = final_data_Spain), na.rm=T)

########################################################################
#                                                                      #
#    1.5 Regression analysis for study1 (aggregated DV)                #
#                                                                      #
########################################################################


## scenario1
summary(lm((agg_s1) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((agg_s1) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)

## scenario2
summary(lm((agg_s2) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((agg_s2) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)

## scenario3
summary(lm((agg_s3) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((agg_s3) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)

## scenario4
summary(lm((agg_s4) ~ (AI.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)
summary(lm((agg_s4) ~ (Human.dummy) + (AI_Human.dummy) , data = final_data_Spain), na.rm=T)


########################################################################
#                                                                      #
#         1.6 Internal consistency check and factor analysis            #
#                             Appendix-5 and 7                         #
########################################################################
# calculating Cronbach's alpha for internal consistency for Spain sample
# use the subset of the questions whose internal consistency is to be calculated

library(psych)

alpha_s1s1 <- final_data_Spain[, 47:50]
alpha(alpha_s1s1)

alpha_s1s2 <- final_data_Spain[, 51:54]
alpha(alpha_s1s2)

alpha_s1s3 <- final_data_Spain[, 55:58]
alpha(alpha_s1s3)

alpha_s1s4 <- final_data_Spain[, 59:62]
alpha(alpha_s1s4)

## internal reliability for study2

alpha_s2s1 <- final_data_Spain[, 63:66]
alpha(alpha_s2s1)

alpha_s2s2 <- final_data_Spain[, 67:70]
alpha(alpha_s2s2)

alpha_s2s3 <- final_data_Spain[, 71:74]
alpha(alpha_s2s3)

# 2.Choose the "right" number of factors to retain
##parallel analysis
library(nFactors)
library(psych)
library(GPArotation)

correl = cor(alpha_s1s1, use = "pairwise.complete.obs")
correl
symnum(correl)

ev <- eigen(cor(alpha_s1s1)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s1),var=ncol(x = alpha_s1s1), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen
nofactors1 <- fa.parallel(alpha_s1s1, fm='ml', fa= 'fa')
sum(nofactors1$fa.values >.7)


s1.fa1 <- fa(r = cor(alpha_s1s1), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa1,labels=names(x = alpha_s1s1),cex=.7, ylim=c(-.1,1)) 
s1.fa1

############### scenario2############
correl = cor(alpha_s1s2, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s1s2)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s2),var=ncol(x = alpha_s1s2), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors2 <- fa.parallel(alpha_s1s2, fm='ml', fa= 'fa')
sum(nofactors2$fa.values >.7)


s1.fa2 <- fa(r = cor(alpha_s1s2), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa2,labels=names(x = alpha_s1s2),cex=.7, ylim=c(-.1,1)) 
s1.fa2

############### scenario3 ############
correl = cor(alpha_s1s3, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis
library(nFactors) 
ev <- eigen(cor(alpha_s1s3)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s3),var=ncol(x = alpha_s1s3), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors3 <- fa.parallel(alpha_s1s3, fm='ml', fa= 'fa')
sum(nofactors3$fa.values >.7)


s1.fa3 <- fa(r = cor(alpha_s1s3), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa3,labels=names(x = alpha_s1s3),cex=.7, ylim=c(-.1,1)) 
s1.fa3


########################## Scenario4 #######################

correl = cor(alpha_s1s4, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s1s4))# get eigenvalues 
ev$values
ap <- parallel(subject=nrow(x = alpha_s1s4),var=ncol(x = alpha_s1s4), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors4 <- fa.parallel(alpha_s1s4, fm='ml', fa= 'fa')
sum(nofactors4$fa.values >.7)


s1.fa4 <- fa(r = cor(alpha_s1s4), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa4,labels=names(x = alpha_s1s4),cex=.7, ylim=c(-.1,1)) 
s1.fa4

#########################Study2- Issue1 ###################

correl = cor(alpha_s2s1, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s2s1)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s2s1),var=ncol(x = alpha_s2s1), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors1 <- fa.parallel(alpha_s2s1, fm='ml', fa= 'fa')
sum(nofactors1$fa.values >.7)


s2.fa1 <- fa(r = cor(alpha_s2s1), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s2.fa1,labels=names(x = alpha_s2s1),cex=.7, ylim=c(-.1,1)) 
s2.fa1

#########################Study2- Issue2 ###################

correl = cor(alpha_s2s2, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s2s2)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s2s2),var=ncol(x = alpha_s2s2), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors2 <- fa.parallel(alpha_s2s2, fm='ml', fa= 'fa')
sum(nofactors2$fa.values >.7)


s2.fa2 <- fa(r = cor(alpha_s2s2), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s2.fa2,labels=names(x = alpha_s2s2),cex=.7, ylim=c(-.1,1)) 
s2.fa2

#########################Study2- Issue3 ###################

correl = cor(alpha_s2s3, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s2s3)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s2s3),var=ncol(x = alpha_s2s3), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors3 <- fa.parallel(alpha_s2s3, fm='ml', fa= 'fa')
sum(nofactors3$fa.values >.7)


s2.fa3 <- fa(r = cor(alpha_s2s3), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s2.fa3,labels=names(x = alpha_s2s3),cex=.7, ylim=c(-.1,1)) 
s2.fa3


#############################################################################################
#                                                                                           #
#                   2  Multilevel models for experiment-1                                   #
#                                                                                           #
#                                                                                           #
#############################################################################################
########################################################################
#                                                                      #
#                  2.1 Creating variables of interest                  #
#                                                                      #
########################################################################



### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
mlm_data_es <- final_data_Spain

#assinging ID to each participant
mlm_data_es$ID <- 1:nrow(mlm_data_es)
mlm_data_es$agg_s1 <- (mlm_data_es$TRUST + mlm_data_es$FAIR + mlm_data_es$AGREE + mlm_data_es$BIASED)/4
mlm_data_es$agg_s2 <- (mlm_data_es$CIVIL_TRUST + mlm_data_es$CIVIL_FAIR + mlm_data_es$CIVIL_AGREE + mlm_data_es$CIVIL_BIASED)/4
mlm_data_es$agg_s3 <- (mlm_data_es$FACTS_TRUST + mlm_data_es$FACTS_FAIR + mlm_data_es$FACTS_AGREE + mlm_data_es$FACTS_BIASED)/4
mlm_data_es$agg_s4 <- (mlm_data_es$DEC_TRUST + mlm_data_es$DEC_FAIR + mlm_data_es$DEC_AGREE + mlm_data_es$DEC_BIASED)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for study1.

mlm_data_es_L <- melt(mlm_data_es, id.vars = c("ID"), measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))
mlm_data_es_L2 <- melt(mlm_data_es, id.vars = c("ID", "condition") , measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))

library(car)
mlm_data_es_L2$scenario <- recode(mlm_data_es_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;
"agg_s3" = 3; "agg_s4" = 4;')

#Create dummy variables for conditions

mlm_data_es_L2$AI <- ifelse(mlm_data_es_L2$condition == 1, 1, 0)
mlm_data_es_L2$AI_Human <- ifelse(mlm_data_es_L2$condition == 3, 1, 0)
mlm_data_es_L2$Human <- ifelse(mlm_data_es_L2$condition == 2, 1, 0)
mlm_data_es_L2$moderate <- ifelse(mlm_data_es_L2$scenario == 1, 1, 0)
mlm_data_es_L2$civil_reminder <- ifelse(mlm_data_es_L2$scenario == 2, 1, 0)
mlm_data_es_L2$present_facts <- ifelse(mlm_data_es_L2$scenario == 3, 1, 0)
mlm_data_es_L2$news <- ifelse(mlm_data_es_L2$scenario == 4, 1, 0)

########################################################################
#                                                                      #
#              2.2 Multilevel models and main analysis                 #
#                                                                      #
########################################################################


###Multilevel Models for testing H1 and H2

##Human and news recommendation as the reference categories (H1, Appendix 8. Table 1) 

model_e1_es <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                     AI*moderate + AI*civil_reminder + AI*present_facts + 
                     AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts,
                   random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e1_es)
anova(model_e1_es)
r.squaredGLMM(model_e1_es)

##Estimated marginal means (Appendix 8. Table 2)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_spain_L2$scenario_rec <- ifelse(mlm_data_spain_L2$scenario == 4, 1, 
                                         ifelse(mlm_data_spain_L2$scenario == 1, 4,
                                                ifelse(mlm_data_spain_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_spain_L2$condition_rec <- ifelse(mlm_data_spain_L2$condition == 2, 0, 
                                          ifelse(mlm_data_spain_L2$condition == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_spain <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                      random =~ 1|ID, data=mlm_data_spain_L2, control = lmeControl(opt = "optim"))
summary( model_e1_spain)
anova( model_e1_spain)
vif.lme( model_e1_spain)

### Getting the estimated marginal means
emm2 = emmeans::emmeans( model_e1_spain, ~ factor(condition_rec) * factor(scenario_rec))
emm2


#AI as the reference category -- RQ1 (Figure 2)

model_e1_es_AI_facts <- lme(value ~ Human + AI_Human + moderate + civil_reminder + news +
                              Human*moderate + Human*civil_reminder + Human*news + 
                              AI_Human*moderate + AI_Human*civil_reminder + AI_Human*news,
                            random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e1_es_AI)
anova(model_e1_es_AI)

## Reestimated models with present facts as the reference category for Figure 2

###Multilevel Models for Figure 2

#Human as the reference category 

model_e1_es_facts <- lme(value ~ AI + AI_Human + moderate + civil_reminder + news +
                           AI*moderate + AI*civil_reminder + AI*news + 
                           AI_Human*moderate + AI_Human*civil_reminder + AI_Human*news,
                         random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e1_es_facts)
anova(model_e1_es_facts)

# AI as reference

model_e1_es_AI_facts <- lme(value ~ Human + AI_Human + moderate + civil_reminder + news +
                              Human*moderate + Human*civil_reminder + Human*news + 
                              AI_Human*moderate + AI_Human*civil_reminder + AI_Human*news,
                            random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e1_es_AI_facts)
anova(model_e1_es_AI_facts)


########################################################################
#                                                                      #
#       2.3   Multilevel models for pre-registered outcomes            #
#                     (Appendix 10. Experiment 1)                      #
#                                                                      #
########################################################################


####################################### FAIR ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1/

mlm_data_spain_fair <- melt(mlm_data_spain, id.vars = c("ID", "condition", "EDU") , measure.vars = c("FAIR", "CIVIL_FAIR", "FACTS_FAIR", "DEC_FAIR"))

library(car)
mlm_data_spain_fair$scenario <- recode(mlm_data_spain_fair$variable, '"FAIR" = 1; "CIVIL_FAIR" = 2;
"FACTS_FAIR" = 3; "DEC_FAIR" = 4;')

#Create dummy variables for conditions

mlm_data_spain_fair$AI <- ifelse(mlm_data_spain_fair$condition == 1, 1, 0)
mlm_data_spain_fair$AI_Human <- ifelse(mlm_data_spain_fair$condition == 3, 1, 0)
mlm_data_spain_fair$Human <- ifelse(mlm_data_spain_fair$condition == 2, 1, 0)
mlm_data_spain_fair$moderate <- ifelse(mlm_data_spain_fair$scenario == 1, 1, 0)
mlm_data_spain_fair$civil_reminder <- ifelse(mlm_data_spain_fair$scenario == 2, 1, 0)
mlm_data_spain_fair$present_facts <- ifelse(mlm_data_spain_fair$scenario == 3, 1, 0)
mlm_data_spain_fair$news <- ifelse(mlm_data_spain_fair$scenario == 4, 1, 0)

#Multilevel Models for figures 
# Human as reference

model_e1_es_fair <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                          AI*moderate + AI*civil_reminder + AI*present_facts + 
                          AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts,
                        random =~ 1|ID, data=mlm_data_spain_fair, control = lmeControl(opt = "optim"))
summary(model_e1_es_fair)
anova(model_e1_es_fair)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_es_fair$scenario_rec <- ifelse(mlm_data_es_fair$scenario == 4, 1, 
                                        ifelse(mlm_data_es_fair$scenario == 1, 4,
                                               ifelse(mlm_data_es_fair$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_es_fair$condition_rec <- ifelse(mlm_data_es_fair$condition == 2, 0, 
                                         ifelse(mlm_data_es_fair$condition == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_es_mfair <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                         random =~ 1|ID, data=mlm_data_es_fair, control = lmeControl(opt = "optim"))
summary(model_e1_es_mfair)
anova(model_e1_es_mfair)
vif.lme(model_e1_es_mfair)

### Getting the estimated marginal means
emm1_fair = emmeans::emmeans(model_e1_es_mfair, ~ factor(condition_rec) * factor(scenario_rec))
emm1_fair

####################################### BIAS ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1

mlm_data_spain_bias <- melt(mlm_data_spain, id.vars = c("ID", "condition", "EDU") , measure.vars = c("BIASED", "CIVIL_BIASED", "FACTS_BIASED", "DEC_BIASED"))

mlm_data_spain_bias$scenario <- recode(mlm_data_spain_bias$variable, '"BIASED" = 1; "CIVIL_BIASED" = 2;
"FACTS_BIASED" = 3; "DEC_BIASED" = 4;')

#Create dummy variables for conditions

mlm_data_spain_bias$AI <- ifelse(mlm_data_spain_bias$condition == 1, 1, 0)
mlm_data_spain_bias$AI_Human <- ifelse(mlm_data_spain_bias$condition == 3, 1, 0)
mlm_data_spain_bias$Human <- ifelse(mlm_data_spain_bias$condition == 2, 1, 0)
mlm_data_spain_bias$moderate <- ifelse(mlm_data_spain_bias$scenario == 1, 1, 0)
mlm_data_spain_bias$civil_reminder <- ifelse(mlm_data_spain_bias$scenario == 2, 1, 0)
mlm_data_spain_bias$present_facts <- ifelse(mlm_data_spain_bias$scenario == 3, 1, 0)
mlm_data_spain_bias$news <- ifelse(mlm_data_spain_bias$scenario == 4, 1, 0)

#Multilevel Models for figures 
# Human as reference

model_e1_es_bias <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                          AI*moderate + AI*civil_reminder + AI*present_facts + 
                          AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts,
                        random =~ 1|ID, data=mlm_data_spain_bias, control = lmeControl(opt = "optim"))
summary(model_e1_es_bias)
anova(model_e1_es_bias)

#Recoding dummy variables for conditions
mlm_data_es_bias$scenario_rec <- ifelse(mlm_data_es_bias$scenario == 4, 1, 
                                        ifelse(mlm_data_es_bias$scenario == 1, 4,
                                               ifelse(mlm_data_es_bias$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_es_bias$condition_rec <- ifelse(mlm_data_es_bias$condition == 2, 0, 
                                         ifelse(mlm_data_es_bias$condition == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_es_mbias <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                         random =~ 1|ID, data=mlm_data_es_bias, control = lmeControl(opt = "optim"))
summary(model_e1_es_mbias)
anova(model_e1_es_mbias)
vif.lme(model_e1_es_mbias)

### Getting the estimated marginal means
emm1_bias = emmeans::emmeans(model_e1_es_mbias, ~ factor(condition_rec) * factor(scenario_rec))
emm1_bias

####################################### TRUST ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1

mlm_data_spain_trust <- melt(mlm_data_spain, id.vars = c("ID", "condition", "EDU") , measure.vars = c("TRUST", "CIVIL_TRUST", "FACTS_TRUST", "DEC_TRUST"))

mlm_data_spain_trust$scenario <- recode(mlm_data_spain_trust$variable, '"TRUST" = 1; "CIVIL_TRUST" = 2;
"FACTS_TRUST" = 3; "DEC_TRUST" = 4;')

#Create dummy variables for conditions

mlm_data_spain_trust$AI <- ifelse(mlm_data_spain_trust$condition == 1, 1, 0)
mlm_data_spain_trust$AI_Human <- ifelse(mlm_data_spain_trust$condition == 3, 1, 0)
mlm_data_spain_trust$Human <- ifelse(mlm_data_spain_trust$condition == 2, 1, 0)
mlm_data_spain_trust$moderate <- ifelse(mlm_data_spain_trust$scenario == 1, 1, 0)
mlm_data_spain_trust$civil_reminder <- ifelse(mlm_data_spain_trust$scenario == 2, 1, 0)
mlm_data_spain_trust$present_facts <- ifelse(mlm_data_spain_trust$scenario == 3, 1, 0)
mlm_data_spain_trust$news <- ifelse(mlm_data_spain_trust$scenario == 4, 1, 0)

#Multilevel Models for figures 
# Human as reference

model_e1_es_trust <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                           AI*moderate + AI*civil_reminder + AI*present_facts + 
                           AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts,
                         random =~ 1|ID, data=mlm_data_spain_trust, control = lmeControl(opt = "optim"))
summary(model_e1_es_trust)
anova(model_e1_es_trust)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_es_trust$scenario_rec <- ifelse(mlm_data_es_trust$scenario == 4, 1, 
                                         ifelse(mlm_data_es_trust$scenario == 1, 4,
                                                ifelse(mlm_data_es_trust$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_es_trust$condition_rec <- ifelse(mlm_data_es_trust$condition == 2, 0, 
                                          ifelse(mlm_data_es_trust$condition == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_es_mtrust <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                          random =~ 1|ID, data=mlm_data_es_trust, control = lmeControl(opt = "optim"))
summary(model_e1_es_mtrust)
anova(model_e1_es_mtrust)
vif.lme(model_e1_es_mtrust)

### Getting the estimated marginal means
emm1_trust = emmeans::emmeans(model_e1_es_mtrust, ~ factor(condition_rec) * factor(scenario_rec))
emm1_trust

####################################### AGREE ################################################################################################
mlm_data_spain_agree <- melt(mlm_data_spain, id.vars = c("ID", "condition", "EDU") , measure.vars = c("AGREE", "CIVIL_AGREE", "FACTS_AGREE", "DEC_AGREE"))

mlm_data_spain_agree$scenario <- recode(mlm_data_spain_agree$variable, '"AGREE" = 1; "CIVIL_AGREE" = 2;
"FACTS_AGREE" = 3; "DEC_AGREE" = 4;')

mlm_data_spain_agree$AI <- ifelse(mlm_data_spain_agree$condition == 1, 1, 0)
mlm_data_spain_agree$AI_Human <- ifelse(mlm_data_spain_agree$condition == 3, 1, 0)
mlm_data_spain_agree$Human <- ifelse(mlm_data_spain_agree$condition == 2, 1, 0)
mlm_data_spain_agree$moderate <- ifelse(mlm_data_spain_agree$scenario == 1, 1, 0)
mlm_data_spain_agree$civil_reminder <- ifelse(mlm_data_spain_agree$scenario == 2, 1, 0)
mlm_data_spain_agree$present_facts <- ifelse(mlm_data_spain_agree$scenario == 3, 1, 0)
mlm_data_spain_agree$news <- ifelse(mlm_data_spain_agree$scenario == 4, 1, 0)

#Multilevel Models for figures 
# Human as reference

model_e1_es_agree <- lme(value ~ AI + AI_Human + moderate + civil_reminder + present_facts +
                           AI*moderate + AI*civil_reminder + AI*present_facts + 
                           AI_Human*moderate + AI_Human*civil_reminder + AI_Human*present_facts,
                         random =~ 1|ID, data=mlm_data_spain_agree, control = lmeControl(opt = "optim"))
summary(model_e1_es_agree)
anova(model_e1_es_agree)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_es_agree$scenario_rec <- ifelse(mlm_data_es_agree$scenario == 4, 1, 
                                         ifelse(mlm_data_es_agree$scenario == 1, 4,
                                                ifelse(mlm_data_es_agree$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_es_agree$condition_rec <- ifelse(mlm_data_es_agree$condition == 2, 0, 
                                          ifelse(mlm_data_es_agree$condition == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_es_magree <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                          random =~ 1|ID, data=mlm_data_es_agree, control = lmeControl(opt = "optim"))
summary(model_e1_es_magree)
anova(model_e1_es_magree)
vif.lme(model_e1_es_magree)

### Getting the estimated marginal means
emm1_agree = emmeans::emmeans(model_e1_es_magree, ~ factor(condition_rec) * factor(scenario_rec))
emm1_agree

#############################################################################################
#                                                                                           #
#                   3  Multilevel models for experiment-2                                   #
#                                                                                           #
#                                                                                           #
#############################################################################################
########################################################################
#                                                                      #
#          3.1 Creating variables of interest                          #
#                                                                      #
########################################################################

### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.
mlm_data_es <- final_data_Spain
### Computing moderators
mlm_data_es$avgppid <- scale((mlm_data_es$PPID_1 + mlm_data_es$PPID_2 + mlm_data_es$PPID_3 +mlm_data_es$PPID_4)/4)
mlm_data_es$ideostr <- scale(abs(mlm_data_es$IDEO_1-5))
mlm_data_es$catalunya_str = scale(abs(mlm_data_es$CATAL_3_1 + mlm_data_es$CATAL_3_2 + mlm_data_es$CATAL_3_3)/3 - 6.33)
mlm_data_es$immig_str = scale(abs(mlm_data_es$IMM_3_1 + mlm_data_es$IMM_3_2 + mlm_data_es$IMM_3_3)/3 - 6.33)
mlm_data_es$unempl_str = scale(abs(mlm_data_es$UNEM_3_1 + mlm_data_es$UNEM_3_2 +  mlm_data_es$UNEM_3_3)/3 - 6.33)
mlm_data_es$counter_pro <- ifelse(mlm_data_es$pro ==1,0,1)

mlm_data_es$issue1str <- mlm_data_es$immig_str
mlm_data_es$issue2str <- mlm_data_es$catalunya_str
mlm_data_es$issue3str <- mlm_data_es$unempl_str

#creating dataset backup
#assigning ID to each participant
mlm_data_es$ID <- 1:nrow(mlm_data_es)
mlm_data_es$agg_s1 <- (mlm_data_es$IMM_CRED + mlm_data_es$IMM_QUALITY + mlm_data_es$IMM_AGREE + mlm_data_es$IMM_BIAS)/4
mlm_data_es$agg_s2 <- (mlm_data_es$BEN_CRED + mlm_data_es$BEN_QUALITY + mlm_data_es$BEN_AGREE + mlm_data_es$BEN_BIAS)/4
mlm_data_es$agg_s3 <- (mlm_data_es$TE_CRED + mlm_data_es$TE_QUALITY + mlm_data_es$TE_AGREE + mlm_data_es$TE_BIAS)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have three issues we
###       consider each issue as level. hence the data has three levels for study 2.

mlm_data_es_L <- melt(mlm_data_es, id.vars = c("ID", "counter_pro",
                                               "avgppid", "ideostr", "catalunya_str", "immig_str", "unempl_str"),
                      measure.vars = c("agg_s1", "agg_s2", "agg_s3"))
mlm_data_es_L2 <- melt(mlm_data_es, id.vars = c("ID", "condition", "counter_pro", "avgppid", "ideostr", "catalunya_str", "immig_str", "unempl_str", "issue1str", "issue2str", "issue3str") 
                       , measure.vars = c("agg_s1", "agg_s2", "agg_s3"))

library(car)
mlm_data_es_L2$scenario <- recode(mlm_data_es_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;"agg_s3" = 3;')

#Create dummy variables for conditions

mlm_data_es_L2$AI <- ifelse(mlm_data_es_L2$condition == 1, 1, 0)
mlm_data_es_L2$Human <- ifelse(mlm_data_es_L2$condition == 2, 1, 0)
mlm_data_es_L2$AI_Human <- ifelse(mlm_data_es_L2$condition == 3, 1, 0)
mlm_data_es_L2$unemployment <- ifelse(mlm_data_es_L2$scenario == 1, 1, 0)
mlm_data_es_L2$immigration <- ifelse(mlm_data_es_L2$scenario == 2, 1, 0)
mlm_data_es_L2$catalunya <- ifelse(mlm_data_es_L2$scenario == 3, 1, 0)

# for automated table creating same name for the variables
# issue3: unemployment, issue1: immigration,  issue2: catalunya
mlm_data_es_L2$issue3 <- ifelse(mlm_data_es_L2$scenario == 1, 1, 0)
mlm_data_es_L2$issue1 <- ifelse(mlm_data_es_L2$scenario == 2, 1, 0)
mlm_data_es_L2$issue2 <- ifelse(mlm_data_es_L2$scenario == 3, 1, 0)


########################################################################
#                                                                      #
#       3.2 MLM models main analysis (Dr. Joao)                        #
#                                                                      #
########################################################################
#Multilevel baseline Model H3 -- is counter seen more negat than pro (Appendix 9 table 1)

model_e2_es_H3 <- lme(value ~  counter_pro + issue2 + issue1 +
                        issue2*counter_pro + issue1*counter_pro,
                      random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e2_es_H3)
anova(model_e2_es_H3)
vif.lme(model_e2_es_H3)


#### Hypothesis testing for figures: retest with catalunya as the ref category.

# Model: M2 (Model only with counter attitudnal participants) (Appendix 9 table 2)
model_e2_es_H4 <- lme(value ~ Human + AI_Human + issue1 + issue2 +
                        Human*issue1 + Human*issue2 +
                        AI_Human*issue1 + AI_Human*issue2,
                      random =~ 1|ID, data=mlm_data_es_L2[mlm_data_es_L2$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_es_H4)
anova(model_e2_es_H4)
vif.lme(model_e2_es_H4)


# Model: M3 (Model M2 with Mix as reference category): retest with catalunya as the ref
model_e2_es <- lme(value ~ AI + Human + immigration + unemployment +
                     AI*immigration + AI*unemployment +
                     Human*immigration + Human*unemployment,
                   random =~ 1|ID, data=mlm_data_es_L2[mlm_data_es_L2$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_es)
anova(model_e2_es)
vif.lme(model_e2_es)

## EXPLORATORY --  is pro/con seen more posit from human/ai/mixed 

model_e2_es_test <- lme(value ~ Human + AI_Human + counter_pro +
                          Human*counter_pro + AI_Human*counter_pro,
                        random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e2_es_test)
anova(model_e2_es_test)
vif.lme(model_e2_es_test)

model_e2_es_test_AI <- lme(value ~ AI + AI_Human + counter_pro +
                             AI*counter_pro + AI_Human*counter_pro,
                           random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))
summary(model_e2_es_test_AI)
anova(model_e2_es_test_AI)
vif.lme(model_e2_es_test_AI)



###  H 7 heterogeneous effects (Appendix 12. Table 1)
### for pol identity strength for figures

model_e2_es_avgppid <- lme(value ~ Human + AI_Human + immigration + catalunya + 
                             Human*immigration + Human*catalunya +
                             AI_Human*immigration + AI_Human*catalunya +
                             avgppid + avgppid*Human + avgppid*AI_Human + 
                             avgppid*Human*catalunya + avgppid*AI_Human*catalunya + avgppid*Human*immigration + avgppid*AI_Human*immigration,
                           random =~ 1|ID, data=mlm_data_es_L2[mlm_data_es_L2$counter_pro == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_es_avgppid)
anova(model_e2_es_avgppid)
vif.lme(model_e2_es_avgppid)

## Ideology strength as a moderator (Appendix 13 table 1)
model_e2_es_ideostr <- lme(value ~ Human + AI_Human + immigration + catalunya + 
                             Human*immigration + Human*catalunya +
                             AI_Human*immigration + AI_Human*catalunya +
                             ideostr + ideostr*Human + ideostr*AI_Human +
                             ideostr*Human*catalunya + ideostr*AI_Human*catalunya + ideostr*Human*immigration + ideostr*AI_Human*immigration,
                           random =~ 1|ID, data=mlm_data_es_L2[mlm_data_es_L2$counter_pro == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_es_ideostr)
anova(model_e2_es_ideostr)
vif.lme(model_e2_es_ideostr)

# H7 Heterogeneous effects issue attitudes (Interacting scenarios with issue att_strt, only for counter condition)
#(Appendix 13 table 2)
model_e2_es_att <- lme(value ~ Human + AI_Human + immigration + catalunya + 
                         Human*immigration + Human*catalunya +  AI_Human*immigration + AI_Human*catalunya +
                         catalunya_str + immig_str + unempl_str + Human*unempl_str + AI_Human*unempl_str +
                         Human*immigration*immig_str + AI_Human*immigration*immig_str + Human*catalunya*catalunya_str + AI_Human*catalunya*catalunya_str,
                       random =~ 1|ID, data=mlm_data_es_L2[mlm_data_es_L2$counter_pro == 1,], 
                       control = lmeControl(opt = "optim"))
summary(model_e2_es_att)
anova(model_e2_es_att)
vif.lme(model_e2_es_att)
########################################################################
#                                                                      #
#    3.3   MLM models per DV for the experiment2 (Appendix 11 table 2) #
#                                                                      #
########################################################################

# creating table per DV for experiment 2 for Spain
mlm_data_es_e2_agree <- melt(mlm_data_es, id.vars = c("ID", "condition", "counter_pro") , measure.vars = c("IMM_AGREE", "BEN_AGREE", "TE_AGREE"))
mlm_data_es_e2_agree$scenario <- recode(mlm_data_es_e2_agree$variable, '"IMM_AGREE" = 1; "BEN_AGREE" = 2; "TE_AGREE" = 3;')

mlm_data_es_e2_agree$AI <- ifelse(mlm_data_es_e2_agree$condition == 1, 1, 0)
mlm_data_es_e2_agree$Human <- ifelse(mlm_data_es_e2_agree$condition == 2, 1, 0)
mlm_data_es_e2_agree$AI_Human <- ifelse(mlm_data_es_e2_agree$condition == 3, 1, 0)
mlm_data_es_e2_agree$unemployment <- ifelse(mlm_data_es_e2_agree$scenario == 1, 1, 0)
mlm_data_es_e2_agree$immigration <- ifelse(mlm_data_es_e2_agree$scenario == 2, 1, 0)
mlm_data_es_e2_agree$catal <- ifelse(mlm_data_es_e2_agree$scenario == 3, 1, 0)

model_e2_es_H4_agree <- lme(value ~ Human + AI_Human + immigration + catal +
                              Human*immigration + Human*catal +
                              AI_Human*immigration + AI_Human*catal,
                            random =~1|ID, data=mlm_data_es_e2_agree[mlm_data_es_e2_agree$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_es_H4_agree)

mlm_data_es_e2_quality <- melt(mlm_data_es, id.vars = c("ID", "condition", "counter_pro") , measure.vars = c("IMM_QUALITY", "BEN_QUALITY", "TE_QUALITY"))
mlm_data_es_e2_quality$scenario <- recode(mlm_data_es_e2_quality$variable, '"IMM_QUALITY" = 1; "BEN_QUALITY" = 2; "TE_QUALITY" = 3;')

mlm_data_es_e2_quality$AI <- ifelse(mlm_data_es_e2_quality$condition == 1, 1, 0)
mlm_data_es_e2_quality$Human <- ifelse(mlm_data_es_e2_quality$condition == 2, 1, 0)
mlm_data_es_e2_quality$AI_Human <- ifelse(mlm_data_es_e2_quality$condition == 3, 1, 0)
mlm_data_es_e2_quality$unemployment <- ifelse(mlm_data_es_e2_quality$scenario == 1, 1, 0)
mlm_data_es_e2_quality$immigration <- ifelse(mlm_data_es_e2_quality$scenario == 2, 1, 0)
mlm_data_es_e2_quality$catal <- ifelse(mlm_data_es_e2_quality$scenario == 3, 1, 0)

model_e2_es_H4_quality <- lme(value ~ Human + AI_Human + immigration + catal +
                                Human*immigration + Human*catal +
                                AI_Human*immigration + AI_Human*catal,
                              random =~1|ID, data=mlm_data_es_e2_quality[mlm_data_es_e2_quality$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_es_H4_quality)


mlm_data_es_e2_cred <- melt(mlm_data_es, id.vars = c("ID", "condition",  "counter_pro") , measure.vars = c("IMM_CRED", "BEN_CRED", "TE_CRED"))
mlm_data_es_e2_cred$scenario <- recode(mlm_data_es_e2_cred$variable, '"IMM_CRED" = 1; "BEN_CRED" = 2; "TE_CRED" = 3;')

mlm_data_es_e2_cred$AI <- ifelse(mlm_data_es_e2_cred$condition == 1, 1, 0)
mlm_data_es_e2_cred$Human <- ifelse(mlm_data_es_e2_cred$condition == 2, 1, 0)
mlm_data_es_e2_cred$AI_Human <- ifelse(mlm_data_es_e2_cred$condition == 3, 1, 0)
mlm_data_es_e2_cred$unemployment <- ifelse(mlm_data_es_e2_cred$scenario == 1, 1, 0)
mlm_data_es_e2_cred$immigration <- ifelse(mlm_data_es_e2_cred$scenario == 2, 1, 0)
mlm_data_es_e2_cred$catal <- ifelse(mlm_data_es_e2_cred$scenario == 3, 1, 0)

model_e2_es_H4_cred <- lme(value ~ Human + AI_Human + immigration + catal +
                             Human*immigration + Human*catal +
                             AI_Human*immigration + AI_Human*catal,
                           random =~1|ID, data=mlm_data_es_e2_cred[mlm_data_es_e2_cred$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_es_H4_cred)




mlm_data_es_e2_bias <- melt(mlm_data_es, id.vars = c("ID", "condition", "counter_pro") , measure.vars = c("IMM_BIASED", "BEN_BIASED", "TE_BIASED"))
mlm_data_es_e2_bias$scenario <- recode(mlm_data_es_e2_bias$variable, '"IMM_BIASED" = 1; "BEN_BIASED" = 2; "TE_BIASED" = 3;')

mlm_data_es_e2_bias$AI <- ifelse(mlm_data_es_e2_bias$condition == 1, 1, 0)
mlm_data_es_e2_bias$Human <- ifelse(mlm_data_es_e2_bias$condition == 2, 1, 0)
mlm_data_es_e2_bias$AI_Human <- ifelse(mlm_data_es_e2_bias$condition == 3, 1, 0)
mlm_data_es_e2_bias$unemployment <- ifelse(mlm_data_es_e2_bias$scenario == 1, 1, 0)
mlm_data_es_e2_bias$immigration <- ifelse(mlm_data_es_e2_bias$scenario == 2, 1, 0)
mlm_data_es_e2_bias$catal <- ifelse(mlm_data_es_e2_bias$scenario == 3, 1, 0)
model_e2_es_H4_bias <- lme(value ~ Human + AI_Human + immigration + catal +
                             Human*immigration + Human*catal +
                             AI_Human*immigration + AI_Human*catal,
                           random =~1|ID, data=mlm_data_es_e2_bias[mlm_data_es_e2_bias$counter_pro == 1,], control = lmeControl(opt = "optim"))
summary(model_e2_es_H4_bias)


########################################################################
#                                                                      #
#   3.4 Composite Marginal means for the experiment2                   #
#                                                                      #
########################################################################

#Recoding scenario instead of using dummy variables to make visualization easier (coding BEN as reference cat).
mlm_data_es_L2$scenario_rec <- ifelse(mlm_data_es_L2$scenario == 3, 1, ifelse(mlm_data_es_L2$scenario == 1, 3, 2))


#Model for H3, H4 and H5 using factor variables instead of dummy (results are the same, only variable names are less clear in output)
#(Appendix-9 table 3)
model_e2_es_emm <- lme(value ~ factor(condition) * factor(scenario_rec) * factor(counter_pro),
                       random =~ 1|ID, data=mlm_data_es_L2, control = lmeControl(opt = "optim"))

### Condition 1 = AI, Condition 2 = Human, Condition 3 = Mixed
### Issue 1 = Immigration, Issue 2 = unemployment, Issue 3 = catalan independence
### Counter_pro 0 = pro-attitudinal, Counter_pro 1 = counter-attitudinal

summary(model_e2_es_emm)
anova(model_e2_es_emm)
vif.lme(model_e2_es_emm)
### Getting the estimated marginal means
emm1 = emmeans::emmeans(model_e2_es_emm, ~ factor(condition) * factor(scenario_rec) * factor(counter_pro))
emm1



# H7 Heterogeneous effects
### Party ID (Appendix 12 table 2)
#simple model
model_e2_es_avgppid <- lme(value ~ factor(condition) * factor(scenario_rec) * avgppid,
                           random =~ 1|ID, data=mlm_data_es_L2[mlm_data_es_L2$counter_pro == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_es_avgppid)
anova(model_e2_es_avgppid)
vif.lme(model_e2_es_avgppid)

### Getting the estimated marginal means.
#NOTE: Because we have a continupus moderator, we need the argument cov.reduce = range
emm2 = emmeans::emmeans(model_e2_es_avgppid, ~ factor(condition) * factor(scenario_rec) * avgppid, cov.reduce = range)
emm2